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A neutral ecology model is simulated on an island chain, in which neighbouring islands can 
exchange individuals but only the first island is able to receive immigrants from a metacommunity. 
It is found by several measures that biodiversity decreases along the chain, being highest for the first 
island. Subtle changes in taxon abundance distributions can be detected when islands in the chain 
are compared to diversity-matched single islands. The results potentially apply to human microbial 
diversity, but highlight the difficulty of using static single-site taxon abundance distributions to 
discriminate between dispersal limitation mechanisms. 
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It has recently been observed that human microbial di- 
versity varies systematically between body sites [H, IH , for 
example pliylogenetic diversity is higher for the palm of 
the hand and the sole of the foot, than for the armpit 
and forehead (Fig. S14 in Ref. Q). A high degree of 
intcr-individual variability is also observed, to the point 
where the characterisation of residual skin bacteria has 
been proposed novel forensic tool The latter, 
in particular, supports the notion that stochastic disper- 
sal limitation may play a significant role in determining 
microbial diversity. Stochastic dispersal limitation is a 
signature element of Hubbell's unified neutral model of 
biodiversity and biogcography [3| , and this motivates the 
question of whether neutral models can be applied to hu- 
man microbial biodiversity and biogeography. This is a 
hard problem and I do not claim to have solved it here. 
Rather, the present study is restricted to exploring the 
role of dispersal mechanisms in the context of neutral 
theory, keeping in mind the possible application to the 
human microbiome. 

The merits of neutral models have been debated ex- 
tensively jsj, and it is far from obvious that they should 
apply to human microbiota 0, However Hubbell's 
neutral model has recently been successfully applied to 
predict microbial diversity in tree holes [1, H]. In this 
context it is important to note that it has been argued 
dispersal limitation is the dominant factor determining 
taxon abundancies [lol . with other neutral model 
ideas, such as the zero-sum constraint (single trophic 
level; community saturation) or speciation by point mu- 
tation, playing a lesser role. 

If taxon abundancies are largely determined by 
stochastic dispersal limitation, then a couple of limit- 
ing hypotheses (Fig. [ij present themselves to explain the 
observed variations in human microbial diversity. The 
first is a variable-immigration-rate hypothesis in which 
different body sites are envisaged as being microbial 'is- 
lands' in contact with a microbial metacommunity but 
effectively isolated from each other. Here variation in di- 
versity corresponds to a variable immigration rate from 
the metacommunity. The second hypothesis is an island- 
chain hypothesis in which it is envisaged that migration 
can take place between islands but, in extremis^ it is only 
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FIG. 1: Putative explanations for a variation in human micro- 
bial diversity, based on dispersal-limitation and the theory of 
island biogeography: (a) variable-immigration-rate hypothe- 
sis and (b) island-chain hypothesis. 



the first island (e. g. the hand) that receives immigrants 
from the metacommunity. In this case one expects that 
diversity should decrease as one moves along the chain 
away from the island in contact with the metacommu- 
nity, due to dispersal limitation. This is confirmed in the 
present study. 

Of course these hypotheses represent limiting cases 
and, if dispersal limitation is relevant, reality probably 
lies somewhere in between. A second question therefore 
is whether one can use taxon abundance distributions 
to distinguish between dispersal mechanisms. Unfortu- 
nately, the present study finds that both hypotheses lead 
to rather similar abundance distributions. When this 
is conflated with other factors, such as deviations from 
neutral model dynamics [Hj, it is probably going to be 
difhcult to distinguish between dispersal mechanisms on 
the basis of static single-site measurements of microbial 
diversity. 

The neutral model has been extensively studied for iso- 
lated islands in contact with a metacommunity ^4. 10. 13l- 
IToj , but only for certain cases has it been solved for mul- 
tiple islands, or 'patches', which are able to exchange 
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individuals [20|, |2l|. In particular, the island chain prob- 
lem has not been solved (i. e. where individuals can mi- 
grate between neighbouring islands but immigration is 
restricted to the first island in the chain). The primary 
goal of the present study is to solve this problem. Al- 
though in principle one can approach this analytically, 
the experience of Valladc and Houchmandzadeh [l^l for 
two islands suggests this will be effectively unmanage- 
able. I therefore approach the problem by means of sim- 
ulations. 

Let me start by summarising the mathematical char- 
acterisation of taxon abundance distributions. Suppose 
there are K taxa and Ni individuals in the i-th taxon 
(i = 1 . . . K), in a population of J = J2i=i individu- 
als. The relative abundance of the i-th taxon is defined 
to be oji = Nil J. The taxon abundance distribution 
is characterised by <j)k, the number of taxa containing 
k individuals. Formally 0^ = X^iLi ^k,Ni where 5n,m 
takes the value unity if n = m and is zero otherwise. 
Given the set of Ni one can easily calculate (fik- One has 
K = X)fc=i 't'k and J = Since no taxon can 

contain more individuals than there are in the commu- 
nity as a whole, (f>k = for k > J. Similarly = 1 if 
and only if all the individuals belong to the same taxon 
(the 'monodominated' state), otherwise cf>j = 0. 

In standard neutral model dynamics, population sizes 
remain fixed (saturated) and are specified at the outset, 
whilst the number of taxa and the number of individuals 
per taxon fluctuates. I adopt the notation of Vallade and 
Houchmandzadeh [l^, [l^ and write (• • • ) to indicate a 
quantity averaged over an ensemble of populations under- 
going neutral model dynamics. The information in (0^) is 
conveniently represented by giving the ensemble-average 
probability pioj) that an individual belongs to a taxon of 
relative abundance lo OIlMllll- For a community of a fi- 
nite size, p{uj) is a discrete array or 'comb' of 5-functions, 
even after ensemble-averaging, since uj can only take on 
discrete multiples of 1/J. However as J — > oo, p{uj) be- 
comes a continuous function. One can show that the 
continuum limit is p{i-u) = lim/_i.oo k{(j)k) where k = loJ 

I shall additionally use several ensemble-average mea- 
sures of diversity. T he p rincipal one of these is the Simp- 
son diversity index [23|, defined for a given set of taxon 
abundancies to be = 1 — ^.-^-^ wf. It is related to 
the second moment of the taxon abundance distribution 
by = 1 — J^^ TlikLi ^'^'t'k- From this it can easily be 
shown that in the continuum limit 



{D)=l- jlLUp{Lu)d,0. 



(1) 



The second diversity measure is the ensemble-average 
number of taxa {K) = X]fc=i('/'fc)- The third is the 
ensemble-average monodominance probability {4>j) — as 
explained above (j)j is 1 or according to whether or not 
all the individuals belong to the same taxon. 

As an order parameter, the Simpson index {D) has 
some advantages over {K) and {<t>,j)'- it remains well de- 
fined in the continuum limit J — > oo, there are some 
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FIG. 2: Steady state island diversity (D) from neutral model 
simulations on islands of various sizes with varying immi- 
gration rates, normalised by the metacommunity diversity 
= 9/{e + 1). The agreement with the theory, Eq. (Q, 
is excellent. Error bars are smaller than the symbols. 



particularly simple theoretical expressions for {D) un- 
der neutral model dynamics, and in particular there is 
a prediction (confirmed by simulation) that {D) fac- 
torises into a product of the metacommunity diversity 
index, and an island factor. Also the Simpson diversity 
index generalises naturally to a measure of /3-diversity 
poj . and to time-series data The index satisfies 

< {D) < 1 — 1/A'. There is a mild disadvantage in 
that {D) loses sensitivity to the underlying abundance 
distribution at the limiting values. 

Let me next summarise neutral model dynamics. In 
the metacommunity it is as follows. An individual is se- 
lected at random, and with probability 1 — v is replaced 
with a copy of another individual drawn at random from 
the metacommunity, or with probability ly is replaced by 
an individual belonging to a new taxon. Thus v is the 
speciation rate. For = the metacommunity even- 
tually falls into a monodominated state, in an ecological 
analog of the Matthew principle [1^ . For > the taxon 
abundance distribution is a balance between speciation 
and extinction. 

An explicit expression for the taxon abundance distri- 
bution in a metacommunity of size Jm has been obtained 
by a number of workers [j, [13, [H, [lB4I3 • Results are 
quoted as metacommunity (subscript 'M') steady-state 
ensemble- averages : 
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where 6 = {Jm — l)i^/(l — ly). 
for Jm > 1 and < 1. It 
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One has w J]\ji/ 
can be shown that 
Jm (an identity), and {K)m = 
The continuum limit of 
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Eq. ^ can be obtained using Stirling's approximation. 
One finds p{oj) = 9[\ — It follows that the meta- 

community diversity order parameter in the continuum 
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FIG. 3: Approach to steady state for the central island of a 
chain of length n islands starting from either a completely 
monodominated state (lower curves), or a flat abundance dis- 
tribution (upper curves; flat means the Ni axe equahsed sub- 
ject to Ni = J). The n — 1 case is for a single island. 
For n > 1 the steady state is weakly dependent on n. 
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This result was noted by He and Hu by analogy to a 
similar problem in genetics (26j . 

Neutral model dynamics on an island connected to the 
metacommunity is as follows. An individual is selected 
at random, and with probability 1 — m is replaced with 
a copy of another individual drawn at random from the 
island, or with probability m is replaced by an individual 
drawn at random from the metacommunity. Thus m is 
the immigration rate. Similar to the metacommunity, the 
island community eventually falls into a monodominated 
state if m = 0, whereas for m > a steady-state taxon 
abundance distribution arises as a balance between im- 
migration and extinction. It is often a very good approxi- 
mation to assume that the island dynamics are decoupled 
from the metacommunity dynamics; in other words the 
metacommunity can be taken to have a static abundance 
distribution. This is because the metacommunity abun- 
dance distribution turns over on a time scale of order 1/v 
whereas the island abundance distribution turns over on 
a time scale of order 1/m, and typically <C m. Here 
time scales are quoted in terms of the number of replace- 
ment steps per individual, since this is expected to be 
proportional to the real elapsed time 

Exact results for the island taxon abundance distribu- 
tion were obtained only recently [13, 113 j althou gh p artial 
results were obtained by previous authors jl3l4l5| . The 
result is 
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where {x)n = r(.T + n)/V{x) is the Pochhammer sym- 
bol, ((I) = r(J + l)/(r(fc + l)r(J - fc + l)) is the bi- 



nomial coefficient, and = m(J — 1)/(1 — m) plays a 
role similar to 9 for the metacommunity. For J ^ 1 and 
TO <C 1, one has ji = Jm. Note that Jm does not fea- 
ture in this expression, in other words the island abun- 
dance distribution is insensible to the metacommunity 
size. This point is discussed more thoroughly by Vallade 
and Houchmandzadeh [l^ . Eq. (U) simplifies in the limit 
k = J to give an expression for the island monodomi- 
nance probability. 



(5) 



This depends strongly on all the relevant parameters and 
vanishes asymptotically for J — oo at fixed fi as ~ 
erieVfilnJ)-". The continuum limit of Eq. (g]) is [l3, 
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Inserting Eq. ^ into Eq. ([T]) gives a simple but to my 
knowledge previously unreported result. 
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Remarkably, as alluded to above, the diversity index fac- 
torises into the product of the metacommunity diversity 
index {D)m = 0/{6 + 1) and an island factor yu/ [fi + l). 

Simulation of neutral model dynamics as summarised 
above is straightforward. I make the assumption that 
metacommunity and island dynamics are decoupled (dis- 
cussed in more detail below). Therefore I generate a 
large number (lO'^-lO^) of metacommunity abundance 
distribution samples for given Jm and 0, equilibrating 
for 10 X J'li/0 replacement stegs between samples to en- 



sure statistical independence [ISj] . I use these samples in 
subsequent island and island chain simulations. As a ref- 
erencepoint, I shall use 6 = 10, motivated by Woodcock 
et al. [9|, and Jm ~ 10^, motivated not so much by time 
scale considerations (see later) but by the requirement 
that Jm ^ J ^ 1 [iMl- Except where otherwise stated, 
averages are over 10^ samples. 

I undertook a number of single island simulations to 
build confidence in the simulation and analysis method- 
ologies. I find excellent agreement between these simula- 
tions and the theoretical predictions for the steady-state 
properties (equilibrating for 10 x J/m replacement steps 
between samples). For example Fig. [5] compares theory 
and simulation results for the Simpson diversity index. 

The island chain simulations arc performed similarly to 
the single island simulations. I introduce an immigration 
rate mi (for the first island) and an inter-island migra- 
tion rate m2. This is illustrated in Fig. [TJb) . Specifically, 
the dynamics are as follows. An individual is selected at 
random. If the chosen individual lies on the first island, 
it is replaced with a copy of another individual on the is- 
land with probability 1 — toi — TO2, an immigrant from the 
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FIG. 4: Aspects of diversity on island cliains of lengtii n: (a) 
diversity index (D), (b) number of taxa {K), and (c) mon- 
odominance probability ((jf>j). The daslied lines in (b) and (c) 
are the {D)-matched single island results for n = 9 (see text). 



metacommunity with probability mi, or a migrant from 
the neighbouring island vifith probability 7712. If the cho- 
sen individual lies on an island interior to the chain, it is 
replaced with a copy of another individual on the island 
with probability 1 — 2m2, or with a migrant from one of 
the neighbouring islands (selected at random) with prob- 
ability 2m2. If the chosen individual lies on the terminal 
island, it is replaced with a copy of another individual 
on the island with probability 1 — m2, or with a migrant 
from neighbouring island with probability m2. Migrants 
are copies of individuals chosen at random on neighbour- 
ing islands. 

Fig. 13] shows that the approach to steady state of a 
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FIG. 5: Diversity index {D), normalised by the theoretical 
metacommunity diversity index {D)m = d/{d + 1), along a 
chain of 71 = 9 islands: (a) varying inter-island migration 
rate 7712 only, (b) varying migration and immigration rates 
together at fixed island size, (c) varying migration and immi- 
gration rates inversely with island size, (d) varying metacom- 
munity diversity parameter. 
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FIG. 6: Abundance distributions for selected islands in an 
island chain of length n = 9: (a) probability that a randomly 
selected individual belongs to a taxon of relative abundance 
(jj, and (b) cumulative distribution function of the same. The 
results (solid lines) are compared to the theoretical expecta- 
tions for (Z))-matched single islands calculated from Eq. ([6]) 
(dashed lines). Results are averages over 10^ samples. 



chain of n islands is slowed by a factor « 1/n compared 
to the single island case. I therefore equilibrate the island 
chains against each metacommunity sample for at least 
10 X rP'J/m replacement steps between samples, where m 
is the smaller of mi and m2 . In terms of the number of 
replacement steps per individual, the island chain relax- 
ation time scale is of the order n/m. Clearly for > 10 
the assumed time scale separation between this and the 
metacommunity relaxation time Xjv JmI^^ ~ 10* is 
faltering. Nevertheless the results are still valid provided 
they can be shown to be unaffected by varying Jm since 
it can be assumed that in reality Jm is much larger than 
10^ 0. To test this, I repeated many of the simulations 
with Jj\/ = 5 X 10"* and Jm = 2 x 10^. I found this made 
no difFerencc to the measured island properties, within 
the statistical errors. 

For a chain of length n = 9, I examined in some more 
detail how the islands recover their steady state diversity. 
The picture is a little mixed. The curves can be fitted 
by (D) = A-hSe"*/^, but not very well, indicating there 
is not a clearly dominant relaxation time. Moreover the 
fitted value of r is affected by whether one starts with a 
monodominated state or a uniform state. What this all 



suggests is that there is a spectrum of relaxation modes, 
which are excited differently according to the initialisa- 
tion protocol, and which are subsequently mixed up by 
the non-linear dynamics. A more detailed exploration of 
this is left for future work. 

Representative steady-state results for the island chain 
simulations are shown in Figs. HHHl The first conclusion 
(Fig.U) is that diversity decreases, by whatever measure, 
as one moves away from the island in contact with the 
metacommunity. Fig. [S] shows how island diversity varies 
with immigration and migration rates m\ and TO2, island 
size J, and the value of Q. Increasing the inter-island 
migration rate m.^ (Fig. [5l[a)) has the effect of increas- 
ing the diversity along the island chain apart from the 
first island. Additionally increasing the metacommunity 
immigration rate ra\ (Fig.[5l^b)) leads to increased diver- 
sity along the whole chain. Fig. Eljc) supports the notion 
that the island diversity is governed by the combinations 
Jmi and Jm2 rather than the individual values of J, 
TOi and TO2; in close analogy to the theory for the sin- 
gle island. Similarly Fig. EJd) strongly suggests that the 
Simpson diversity index continues to be factorisablc into 
the metacommunity diversity index multiplied by a con- 
tribution from the structure of the island chain, again in 
close analogy to the single island result. 

I next compare islands in the chain to '(D)-matched' 
single islands. Here (D)-matching means a value for /i 
is inferred from Eq. ^ {i.e. [i = {D)/{{D)m - {D)) 
where {D)m = 9/{6 + 1)), and used to calculate values 
of {K) = T,'LMk) and from Eqs. Q and ©. I 
assume the island size J is fixed. The procedure amounts 
to matching the first and second moments of The 
dashed lines in Fig. IH^b) and (c) show systematically that 
the ensemble-average number of taxa is reduced and the 
monodominance probability is increased, comparing an 
island in the island chain with its (D)-matched single 
island counterpart. Thus there is a tendency towards 
fewer, larger taxa, when islands in a chain are compared 
to (D)-matchcd single islands. 

A more detailed examination of the abundance distri- 
butions shows that there is a subtle and non-trivial redis- 
tribution of the taxon abundancies. When compared to 
the (£')-matched single islands. Fig. [5] shows that p(w) 
is reduced for uj < 0.2 and w > 0.8, but increased for 
0.2 < UJ < 0.8. This means that the number of taxa with 
intermediate abundancies is increased at the expense of 
the very rare taxa and the high abundance taxa. But, 
in addition, the cumulative distribution function jumps 
up at w = 1, as shown clearly in Fig. El^b). This corre- 
sponds to the increased monodominance probability. At 
first sight this is at odds with with the redistribution 
towards mid-range abundancies, nevertheless it is a real 
effect and indeed is the reason why monodomination was 
separately studied. 

The loss of the very rare taxa can perhaps be at- 
tributed to the filtering properties of the island chain. 
These taxa are already rare in the metacommunity and 
it could simply be that a representative from a rare taxon 
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is less likely to arrive via migration along an island chain 
than via direet immigration from the metaeommunity (at 
matched (D)). The loss of the high abundance taxa and 
the increased monodominance probability are more mys- 
terious and I do not at present have a clear mechanistic 
explanation. Possibly what is happening for islands in 
a chain, compared to (D)-matched single islands, is that 
the monodominatcd state {uj = 1) has become 'stickier' in 
dynamical terms, without actually becoming an adsorb- 
ing state. In the monodominatcd state there is of course 
only one taxon, with w = 1, and this may come at the 
expense of the high abundance taxa with 0.8 < cj < 1. 

For (D) < 0.5 the abundance distribution for a (D)- 
matched metaeommunity is almost exactly the same as 
that for a (I?)-matched single island. By this I mean 
that p{uj) = 61(1 -w)"-! with 61 = (D)/(l -(£>)) is a very 
good approximation to p(w) from Eq. ([6]). However a 
complete comparison with an equivalent metaeommunity 



is frustrated by the residual dependence of (K) m and 
{4>j)m on the metaeommunity size Jj\/. 

Despite these subtleties, it is clear from Fig. |6] that 
the taxon abundance distributions on an island chain 
arc quite well approximated by (Z3)-matchcd single is- 
lands. This is the origin of the claim in the introduction 
that it is probably going to be difficult to distinguish be- 
tween dispersal mechanisms on the basis of static single- 
site measurements of the taxon abundancies. To resolve 
this question, or indeed to distinguish between dispersal- 
limitation and niche-adaptation (3] , probably requires 
more detailed examination of the /3-diversity [lO, [U, HI] , 
and dynamics [H, [H, [l^l- In this context it may be 
useful to explore spatial correlations [l^ and temporal 
correlations [2J| , which are natural generalisations of the 
Simpson index D. 

I thank Mike Gates, Chris Quince, Bill Sloan, and Dave 
Taylor, for helpful comments. 
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